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Abstract 



The paper is devoted to the efficient computation of high-order cubature 
formulas for volume potentials obtained within the framework of approximate 
<^ approximations. We combine this approach with modern methods of struc- 

tured tensor product approximations. Instead of performing high-dimensional 
discrete convolutions the cubature of the potentials can be reduced to a certain 
number of one-dimensional convolutions leading to a considerable reduction 
O °f computing resources. We propose one-dimensional integral representions 

of high-order cubature formulas for n-dimensional harmonic and Yukawa po- 
tentials, which allow low rank tensor product approximations. 



1 Introduction 



The construction of efficient representations of multi-variate integral operators plays 
a crucial role in the numerics of higher dimensional problems arising in a wide range 
of modern applications. Let us mention multi-dimensional integral equations and 
volume potentials of elliptic and parabolic partial differential operators in W 1 , n > 3. 

In the present paper we study the combination of high-order semi-analytic cu- 
bature formulas for volume potentials with modern methods of structured tensor 
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product approximations. The cubature formulas have been obtained in 0, [6] using 
the method of approximate approximations, see also [7J and the reference therein. 

The application of tensor product approximations to the approximation of volume 
potentials is described for example in [3J IH 12] ■ The main idea is to derive accurate 
tensor product approximations of the density u and the kernel g of the convolution 
integral of the form 

R n 

p=i j=i 

where x = (xi, . . . , x n ) and 

R n 
q=l j = l 

such that the convolution integral can be approximated by one-dimensional convo- 
lutions 

/R n „ 

g(x - y)u{y) dy^Y b i u p II / 9j q \ x j - V^iv) d V ■ 

Then the numerical computation of the integral does not require to perform an 
n-dimensional discrete convolution, for example, instead one has to compute pq 
one- dimensional discrete convolutions, which can lead to a considerable reduction 
of computing time and memory requirements, and gives the possibility to treat real 
world problems. 

In this note we present some variants for the tensor product approximation of 
high order cubature for harmonic, Yukawa and heat potentials as integral operators 
in M. n . After a brief introduction into cubature formulas based on approximate 
approximations and their error behavior, we describe their representations as one- 
dimensional integrals suitable for tensor product approximation. For the example 
of the harmonic potential in IR n , 3 < n < 6, we report on numerical tests for second 
and fourth order formulas, which provide estimates of the tensor rank required to 
approximate the action of the potential on a basis function with a prescribed relative 
error. Further we report on numerical tests for the Yukawa potential, which can be 
approximated accurately with a very small number of tensors in a wide range of 
arguments. In the final section [5] we study the volume heat potential in M. n , derive 
approximation results and describe its tensor product approximation. 



2 Semi-analytic cubature formulas for potentials 

Here we collect some results on high-order cubature formulas for the volume poten- 
tials of the differential operators —A and —A + a 2 , Re a 2 > 0, in M. n . 
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2.1 Harmonic potentials in M. 71 

The harmonic potential is the inverse of the Laplace operator and given in M. n by 
the formula 

The volume potential provides the unique solution of the Poisson equation 

-Af = u inM n , |/(x)| < C|x| n ~ 2 as |x| -> oo . 

The theory of approximate approximations proposes semi-analytic cubature formulas 
for harmonic potentials by using quasi-interpolation of the density u by functions 
for which the integral operator can be taken analytically. For example, approximate 
u by the quasi-interpolant 



(x)=p-»/ 2 53«(/im)^(^^) ) (2.2 



'x — hm s 

u h W = U—'* 2^ Hhm)rjl—=- 
mez" visa 

with a suitable generating function r\. Then the sum 

h 2 \ /x - /im\ 

^nMh(x) = 2^ j (2.3) 

is a cubature of the harmonic potential. If one wants to compute the harmonic po- 
tential of u on the given grid {hk} then one has to compute the discrete convolution 

h 2 

C n u h {hk) = — — fl k-m u(hm) , keZ n , (2.4) 



£>n/2-l 

meZ" 



with the coefficients 



a k = a k (D) = C nV ( *=) = T{ ~ 2 ,} ] [ 0^ dy . (2.5) 



It has been shown that for sufficiently smooth and compactly supported functions the 
cubature formula (2.3) provide approximations with the error 0(h 2M ) + 0(e~ v * h 2 ) 
if the function 77 is chosen as 

r ?2M (x) = vr-"/ 2 4'; / _ 2 1 ) (|x| 2 )e-l x l 2 , (2.6) 
where L ? are the generalized Laguerre polynomials 

Additionally, there holds the analytic representation 

/k\ 1 D"/ 2 - 1 (n |k| 2 x e-^lv^Lf'^^/V) f 
^ 2M ^j=^^ = , 7 (--l,— + . (2.7) 
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Here 7 is the lower incomplete Gamma function defined by 



7(0,2;) 



T a ~ 1 e- T dr. 



(2i 



l(k,x) = (k - 1)!(1 



In the case k G N 



whereas for odd space dimension n one can use 

1 



fc-i 

£ 



;r-' 



71^ 



erf (y/x) 



with the error function erf and the recurrence relation 



7(n/2 + l,:r) 



7(71/2,2;) -e" x x n/2 



to derive analytic expressions for C n V2M- 

The asymptotic error estimate 0(h 2M ) + 0(e~' D7r2 /i 2 ) for the cubatu re fo rmula 



(2.3) is based on the error estimate 0(h 2M ) + e for the quasi-interpolant (2.2) if the 
sufficiently smooth and decaying generating function 77 is subject to the moment 
condition 



77(x)dx = l, / x"r/(x) rfx = 0, Wa , 1 < |a| < 2M . 



(2.9) 



The saturation error e does not converge to zero for h — > 0, but because of 



e = O ( max | ^ Tr\{ y/Vk) e 



27ri(k,x)//i 



keZ n \{0} 



(2.10) 



can be made arbitrarily small if the parameter T> is sufficiently large . Ad dition- 
ally, the harmonic potential maps the fast oscillating saturation term (2.10) into a 



function with norm of order 0(h 2 e), which establishes the error estimate^ h 2M ) + 
0(e~ Vn h 2 ) for the approximation of the harmonic pote ntial by using quasi-interpo- 
lation of the density with the generating function (2.6). 



So in numerical computations with T> > 3 the formulas (2.4 2.7) behave like high 
order cubature formulas for harmonic potentials in R n . The approximation of the 
potential on the grid {hk.} can be done by fast convolutional methods, but one has 
to store the values {a^} for a quite large subset of k G Z". 



2.2 Yukawa potentials in M. n 

The fundamental solution of the operator —A + a 2 , Re a 2 > in M n is given as 



(27r) n / 2 V a 



X|\W2. 



K, 



ra/2-1 



>l x l) 
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where K v is the modified Bessel function of the second kind, also known as Mac- 
donald function, [H 9.6]. Thus the volume potential 



J « a (x-y)«(y) dy 



for n = 3 also called Yukawa potential, provides the solution of the equation 

{-A + a 2 )f = u inM n . 

To derive a cubature formula for that potential we look for a solution of 

-A/ + a 2 / = e" |x|2 , xGf, 
which is given as the one-dimensional integral 

— |x| n/2—1 f 

/(x) = 6 " / ^ n/2 _i(ar)/ n/2 _ 1 (2|x|r)r e"' 2 dr, (2.11) 



o 



where I;, is the modified Bessel function of the first kind, see [5 Section 5.2]. 

Using the known analytic e xpress ions of I n+ i/2 and K n+1 / 2 (cf. pQ) it is possible 
to derive analytic formula of ( |2.11 ) for odd space dimension n. In particular, if 



n = 3, then 

v^e-W 2 / (,a 



/(*) = X ~ |X|) ) " W ^ + |X|) )) ' (2 ' 12) 

where w denotes the scaled complementary error function 



w(z) = e z2 



z 

erfc(-iz) = e-* 2 (l + ^ ^ / dt), (2.13) 



and 

erfc(r) = 1 - erf(r) (2.14) 
is the complementary error function. 

Using the representation from \7\ Theorem 3.5] 



one can derive as in the case of harmonic potentials semi-analytic cubature formulas 
for the Yukawa potential with the approximation rate 0(h 2M ) + 0(e~ V7T h 2 ). 
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3 Tensor product expansions of potentials acting 
on Gaussians 



To obtain a tensor product approximation of the second order cubature formula 

for the harmonic potential with 

x — hm , / i 1 2 \ . « 1 fn 



and C n (e I ' I 2 ) ( x ) = _^_ 7 Q _ 1 ; | x | 



Vh " 11 

we use the formula obtained in [8] 

I r -|x|V(l+t) 

v n J A J (1 + t)™/2 



which is valid for n > 3. The simple quadrature of the integral 

« e -|x| 2 /(i+^) « " e -^ 2 /(i+^) 
£„(e-H )(x)«E-»-(I T ^ = E^n(I T ^ F i 

k=l k=l j=l 

with certain quadrature weights u>k and nodes gives already a tensor product 
approximation. Hence, for x ) and m = (mi, . . . , m n ) 

'X-/im\ £>/l 2 ™ e -{xj-h mj f/{Vh 2 {l+T k )) 



(l + r fc )"/ 2 



£ " (eH ' |2 K^) K 4p^7 5 5:^n 

v v ' fc=l .7 = 1 

which implies that one can approximate 

Vh 2 A JLe-^-m^Kvii+ru)) 
C n u h (hk) « /2 ^ufam) 2^*11 , 

v ; meZ" fe=i i=i v fe; 

k = (ki, . . . , k n ). 

To obtain a similar tensor product approximation for higher order cubature for- 
mula we note that one obtains the s ame convergence order 0(h 2M ) + 0(e~ V7T h 2 ) as 
in the case of generating functions (2.6) if the density is approximated by the sum 

Mm«(x) = (ttX?)-"/ 2 u(hm)f[ri 2M ( Xj ~ hm > 



where the generating function is the tensor product of the one-dimensional generat- 
ing functions 

V 2M (x)=L^\(x 2 )e~ x2 . (3.2) 



6 



and obviously satisfies the moment condition (2.9). 

To get the one-dimensional integral representation of C n ^ rij=i VzmJ we use the 
relation 



M ~ l I i\k J2k lvi-l , , K m . 



fc=0 fc=0 k=0 

where Hk(x) denotes the Hermite polynomial 



H k (x) = (-l) k e< B -[^ \ e 



dx 



k 

2 



Then the solution of the Poisson equation 



-Au(x) = JJ^mO^ 

3=1 



is given by the integral 



a 11 1 tufc / n , +u/2 dt 



i/n(E 



A;!4 fc dxf J (1 + t) 1 / 2 
' o 

(-l) fe d 2fc e _ x .2 /(1+t) ^ dt 



k\ 4 fc ctef / (1 + t) n / 2 

1 °? n M-l 2 

. | ne -^ (E ^_ 4 -v, ( ^ ))d , 

£ 3=1 k=0 V ; 

Again, a tensor product representation of this integral and consequently of the con- 
volution matrix for the high order cubature of the harmonic potential is given by a 
quadrature of the last integral. 

So the problem is reduced to find efficient quadrature formulas for the parameter 
dependent integrals 

|xf/(l+t) 7 " P- X i/( 1+ *) 



(3.3) 



7 i( x ) = / V "x io dt= IT\ 6 / dt . 



°S n M-l 2 

/ M (x) =/ n e-^" (E IT T^ i " ,/2, (lTi))*- 

3 =1 ^=0 V ^ 



More precisely, one has to find a certain quadrature rule with minimal number of 
summands which approximates the integrals with prescribed error for the parameters 
Xj = (kj — rrij)/\/V within the range \xj\ < K and some given bound K. 
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3.1 Quadratures 



It is well known that classical trapezoidal rule is exponentially converging for certain 
classes of integrands, for example periodic functions and rapidly decaying functions 
on the real line. For example, Poisson's summation formula yields that 

oo oo „ 

h E /(**) = E /(it) ( 3 - 4 ) 

k=— oo j=—oo 

for any sufficiently smooth function, say of the Schwarz class <S(R). Here / is the 
Fourier transform 



/(A)= J f(x)e~ 2 ™ x dx. 

— oo 

Thus, 

oo 

J f( x )dx-h e /(^) = E/(ir)' 

_ 00 fc=-oo jjiO 

which indicates that by choosing special substitutions such that the integrand trans- 
forms to a rapidly decaying function with rapidly decaying Fourier transform, the 
trapezoidal rule of step size h can provide very accurate approximations of the in- 
tegral. 

Here we follow a proposal made by J. Waldvogel [9] to compute accurately inte- 
grals of analytic functions. We make the substitutions 

t = e ? , £ = a(r + e r ) and r = b{u - e -tt ) , (3.5) 



where a, b > are certain constants. Then the integrals (3.3) are transformed 
to integrals over K of doubly exponentially decaying integrands /, i.e. |/(w)| < 
cexp(— aexp(|w|)) for \u\ — > oo with certain constants c, a > 0. It is known (cf. e.g. 
[I]) that by suitable truncation of the infinite sum in the trapezoidal rule of step 
size h provides exponentially convergent numerical quadrature algorithms with the 
error estimate 0(e~ c ^ h ). 

After the substitution we have 

oo 

f / |x| 2 x (1 + e-")(l + exp(b(u - e~ u ))) <j>(u) 
7 1 (x) = a6 / exp - — — — — du, 



l + <f)(u)J (l + 0(n))™/ 2 

— oo 

where we set 

(f>(u) = exp(ab(u — e~") + aexp(b(u — e _u ))) . 

Similarly 

Jm(x) = ab j ^\gj{u)du with the functions 

•i_ .7=1 



#(ti) = (1 + e-)(l + exp(% - e-))) </>(«) exp(- E Q 



i r/-Va) 



2 M-i r. , 



1 + 0(«)/ ^ (1 + 0( M )) fe + 1 /2 • 
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The integrals are approximated by the finite sum 



I f(u,x)du^h f(hk,x), |x| < K. 

J 1 A T . 



(3.6) 



k=-N 



3.2 Numerical Results 



3.2.1 Approximation to the integral Ji(x) 



We assume in (3.5) a = b = 1. Figure IT] illustrates the graph of the integrand 
function f(u, x), w G (—4,4), n = 3, for different values of |x| < 10 3 . A similar 
behavior holds for n = 4, 5, 6. 




Figure 1: The plot of the integrand function f(u, x) (a = b = 1) in ii(x) for 
|x| = 0, 1, 10, 100, 500, 1000 (from the left to the right) in the interval u G (-4, 4). 

Table [T] presents the maximum step h and the minimum number of quadrature 
points required to achieve the relative error e, uniformly in |x| G [0, 10 3 ]. We have 
considered the space dimension n = 3,4,5,6. 

It is possible to play with different parameters a and b in order to diminish the 
number of summands in the quadrature formula. Consider e.g. the case a = 6 and 
b = 5. Figure [2] shows the graph of f(u, |x|), u G (0,0.85) for different values of |x|. 
The numericalresults for this quadrature are given in Table [2] 

3.2.2 Approximation to the integral -?2( x ) 
Next, we discuss the computation of the integral 



n 



e '+ 1 



+ 



y — 



y/t + 1 (t + 1) 3 / 2V 2 t+1 



dt. 
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n 


= 3 


Relative 


h 


at i r 

IN umber ol 


Error 




quadrature points 


10- 1 


0.264 


18 


t n— 3 
1U 


n 1 Q7 
U.lo ( 


do 


10~ 5 


0.072 


82 


i n— 7 
ID 


0.055 


lib 


i n— 9 
1U 


0.04o 


lol 


1 n— 11 
1U 


0.036 


on c 

205 




n 


= 5 


Relative 


i 

h 


AT "U £ 

JN umber ol 


Error 




quadrature points 


KT 1 


0.181 


21 




0.088 


59 


KT 5 


0.065 


83 


io~ 7 


0.060 


96 


1Q -9 


0.037 


169 


10" 11 


0.033 


200 


Table 1 


: The approximation of Ji(> 




n 


= 3 


Relative 


h 


Number of 


Error 




quadrature points 


io- 1 


0.0297 


10 


1 n— 3 


0.0125 


28 


10~ 5 


0.0077 


61 


1 n— 7 
11) 


0.0055 


111 


1 n— 9 
10 


0.0042 


170 


1 n— 11 
1U 


0.0034 


247 




n 


= 5 


Relative 


ho 


Number of 


Error 




quadrature points 


io- 1 


0.0380 


7 


IO" 3 


0.0120 


27 


IO -5 


0.0069 


57 


IO" 7 


0.0046 


112 


IO" 9 


0.0034 


179 


IO' 11 


0.0031 


221 





n 


= 4 


Relative 


h 


AT "U £ 

IN umber ol 


Error 




quadrature points 


io- 1 


0.198 


20 


1 n- 3 
1U 


u.uyo 


5z 


10 


0.072 


77 


i n-7 
1U 


U.Uol 


1 Ol 

Izl 


1 n— 9 
1U 


U.U4U 


lo4 


1 n— 11 
1U 


U.Uoo 






r?, 


= 6 


Relative 


/io 


AT "U -C 

IN umber ol 


Error 




quadrature points 


io- 1 


0.156 


26 


io- 3 


0.090 


55 


IO" 5 


0.059 


90 


10~ 7 


0.044 


130 


10~ 9 


0.035 


178 


io- 11 


0.029 


220 



with a = b = 1 in (3.5). 



n 



Relative 
Error 


h 


Number of 
quadrature points 


io- 1 


0.0234 


10 


IO" 3 


0.0107 


30 


10~ 5 


0.0070 


58 


IO" 7 


0.0049 


107 


10~ 9 


0.0037 


169 


io- 11 


0.0033 


217 


n = 6 


Relative 
Error 


h 


Number of 
quadrature points 


io- 1 


0.0185 


12 


IO" 3 


0.0083 


36 


IO" 5 


0.0058 


70 


IO" 7 


0.0042 


117 


IO" 9 


0.0037 


158 


io- 11 


0.0028 


242 



Table 2: The approximation of ii(x) for |x| < 10 3 , with the choice a = 6; b = 5 in 
(pb. 
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Figure 2: The plot of the integrand function f(u, x) (a = 6, b = 5) in Jj(x) for 
|x| = 0, 1, 10, 100, 500, 1000 (from the left to the right) in the interval u G (0, 0.85). 



using the variable transformations (3.5) and the trapezoidal rule (3.6), for n = 3 
and n — 4. In the numerical results below, for the sake of simplicity, we assumed 
x = (x, x, x), with |x| < 10 3 . 




Figure 3: The plot of the integrand function f(u, x) in ^(x) for |x| = 
0,1,10,100,500,1000 (from the left to the right) in the interval u E (-4,4), 
a = 1, b = 1. 



Numerical results for this quadrature are presented in Table [3j with the parame- 
ters a = b = 1, and in Table [4] in the case a = 6, b = 5. 
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n 



n 



Relative 
Error 


n 


Number of 
quadrature points 


10- 1 


0.295 


16 


10" 3 


0.133 


40 


10" 5 


0.072 


82 


10- 7 


0.055 


118 


10- 9 


0.043 


163 


10- 11 


0.036 


204 


Table 3: The approximation of I 
n = 3 


Relative 
Error 


h 


Number of 
quadrature points 


KT 1 


0.0297 


10 


10~ 3 


0.0125 


30 


10- 5 


0.0077 


63 


10- 7 


0.0052 


114 


10- 9 


0.0042 


175 


10' 11 


0.0034 


234 



-tveiative 


h 


INUHlDei 01 


Error 




quadrature points 


10- 1 


0.198 


23 


10~ 3 


0.095 


52 


10~ 5 


0.072 


77 


10- 7 


0.051 


121 


10~ 9 


0.040 


163 


10- 11 


0.033 


206 



'x) for Ixl < 10 3 , with a = b = 1. 



n 



Relative 


h 


Number of 


Error 




quadrature points 


10- 1 


0.0197 


11 


10' 3 


0.0107 


30 


10- 5 


0.0074 


57 


10- 7 


0.0046 


120 


10- 9 


0.0037 


175 


HT 11 


0.0033 


222 



Table 4: The approximation of -^(x) for |x| < 10 3 , with a = 6; b = 5. 



4 Yukawa potential 

To derive a tensor product approximation of the second and higher order cubature 
formulas for the Yukawa potential we use the relation 

1 7e~ a2 */ 4 e~ |x|2/(1+i) 
« a (x-y)e-W dy=-J (1 + t)w/2 dt, 

R" 

obtained in [8] and is valid for all n > 2, see also [7> Theorem 6.4]. Hence, an 
approximate solution of the equation in M. n 

-Af + a 2 f = u 

is given by 



00 



Dtf r e -a 2 Vh?t/A e -|x-/im| 2 /(D/i 2 (l+t)) 

fhM = -r, ^ /n 7 u(hm) / - dt , 

v 1 4(vr£>) n / 2 ^ n K 1 J (l+t) n / 2 



me 





2^ 



which converges with the order 0(h) to /. 
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Analogously to the case of harmonic potentials we consider the integral 
J «a(x - y) J] ^m(%) dy = Y[[J2 dx^) J Ka(x ~ y) e_ly ' 2 dy 



n( V ( ' — / e" a2 '/ 4 6 
, A;!4 fc cfef 7 (1 + t) 1 / 2 



o 



3=1 fc=0 



which is the basis of the cubature formula of the order 0(h 2M ) + 0(e Vn2 h 2 ). 
To get doubly periodic integrands for the integrals 



r e -a 2 t/4 e -W 2 /(l+t) r " e -^ 2 /(l+t) 

^(x) = / 6 - 6 W9 dt = / e"^ 4 FT 6 . dt 

o o 

oo 



3=1 fc=0 V ' 



n M-l , _ , r 2 

j=l fc=0 

we make the substitutions 

t = exp(6(it — exp(— u)) , b > 
and apply the trapezoidal rule to 



oo 

Ki(x) =a J exp(- - 



|x| 2 c 2 Mu) \ (1 + e~ u ) (f)(u) 



+ cj){u) 4 / (1 + 0(m))™/ 2 



/II 
e"^ w/4 (l + e~ u ) 4>(u) \{ exp (- - 



<j>(u 

(1 + cf){u)) k + l l 2 
4.1 Approximation to the integral ifi(x 



M-l 1 2 

X V 1 r(-l/2)f X J \ ^ 

fc=0 



We apply the quadrature formula (3.6) to the integral Ki(x) for n = 3, in the cases 
a 2 = 0.01, a 2 = 0.1 (Table g, a 2 = 1 and a 2 = 4 (Table g. We assumed 6=1. 
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8x 


10" 






6x 


JO- 






4x 


ID" 






2x 


10" 







-3 -2 -1 



0. 15 ■ 
0. 125 ■ 
0.1- 
0. 075 ■ 

0. 05 ■ 
0. 025 ■ 

-1 

2.5x 10' 
2x 10 

1.3k 10' 
lxlO"" 
5x 1 0-" ' 



-1 



Figure 4: The plot of the integrand function /(w,x) in i£i(x), a 2 = 0.01,6 
|x| = 0, 1, 10, 100, 500, 1000 (from the left to the right) in the interval u G (- 



= 1 for 
3,3). 




4x 10' 
3x 10" 
2x 10" 
lx 10" 



1 2 



Fij 
|x 



jure 5: The plot of the integrand function f(u, x) in i^i(x), a 2 = 4, 
= 0, 1, 10, 50, 100, 200 (from the left to the right) in the interval u E ( 



1 2 

6=1 for 
-2,2). 



5 Heat potential 

Consider the non-homogeneous (linear) heat equation 

f t -vA x f = u(x,t), xGl", t>0 
with the initial condition 

/(x,0) = 0, xel". 
It well known that the solution of this Cauchy problem can be written as 

t 



/(x,t) = yW A u(-,A))(x) dX, 



(5.1) 
(5.2) 

(5.3) 
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= 0.01 



= 0.1 



rvcidti vt. 


IbQ 


NiiYnhpT 1 nf 
IN lililUt;! Ul 


±1/1 1 UI 




L^LldU.1 dlj 111 C JJL»lllLo 


io- 1 


0.99 


9 


IO" 3 


0.72 


15 


10~ 5 


0.58 


20 


io- 7 


0.47 


25 


io- 9 


0.39 


32 


10" 11 


0.30 


43 


io- 13 


0.26 


50 


IO" 15 


0.25 


56 



\-< o 1 o "f" n rn 

rvcldll Vc 


h 


IN UIIlUcI Ul 


±L1 1 (JI 




m_idU.IdLUIc pUlllLo 


lO^ 1 


0.98 


7 


IO" 3 


0.68 


12 


IO" 5 


0.58 


17 


io- 7 


0.42 


16 


io- 9 


0.40 


25 


10" 11 


0.29 


36 


IO" 13 


0.25 


43 


IO" 15 


0.23 


50 



Table 5: The approximation of K 1 (x) for |x| < 10 3 , with a 2 = 0.01 (on the left) and 
a 2 = 0.1 (on the right). 



a 2 = l 



a 2 = 4 



Relative 


h 


Number of 


Error 




quadrature points 


io- 1 


0.99 


6 


io- 3 


0.65 


10 


10~ 5 


0.48 


15 


io- 7 


0.38 


20 


10~ 9 


0.37 


22 


10" 11 


0.29 


28 


IO" 13 


0.25 


34 


io- 15 


0.21 


42 



Relative 


h 


Number of 


Error 




quadrature points 


io- 1 


0.92 


5 


io- 3 


0.58 


9 


10~ 5 


0.44 


13 


io- 7 


0.36 


17 


IO' 9 


0.31 


21 


IO" 11 


0.27 


25 


IO" 13 


0.25 


29 


io- 15 


0.16 


46 



Table 6: The approximation of -^i(x) for |x| < 10 3 , with a 2 = 1 (on the left) and 
a 2 = 4 (on the right). 



where Vt is the Poisson integral 



(5.4) 



An approximation of this solution /(x, t) can be obtained if the function / is 
approximated by the quasi-interpolant on the rectangular grid {h m,rj), with h > 
and r > 0, 



-(n+l)/2 



-(t-Tjf/(V T 2 ) e -\*-h m \ 2 /(Vh 2 ) 



(5.5) 



m£Z" 
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Then the sum 

t 

A,r(x,t) = J{V t -AU h>T {;\)){x.)d\ (5.6) 





-n/2 _ \ -(X- T j) 2 /(V r 2 ) 



Is <hm,Tj) J - _ A)n/2 dX 

x / e -|x-e| a /(M*-A))-|€-ftm|V('» 3 ^)^ 



m£Z" 



ftn ^ e -(A-(t-Tj)) 2 /(D r 2 ) e -[x-/i.m| 2 /(d 2 +^A) 

> u(hm,rj) / — — — r— -rz dX 



tt^ 1 )/ 2 ^ K ' J ' J (Vh 2 + 4uX) n / 2 



provides an approximation of /(x, t). 

Since Uh tT (x,t) approximates w(x, t) with 



\u 



(x,t)-u hiT (x,t)| <£ + c((r v / ^) 2 + (/iV / P) 2 ), VxGr,tG[0,T] (5.7) 



the function /^(x, t) approximates the solution /(x, i) with the error 



1 r /■ e -|x-4l 2 /(4(<-A)) 



|/(x,t)-/ft, T (x,t)| = (4?r ^ ra/2 y dA y ( t -A)"/2 N& A )- M ft,r(g,A)|<% 

R™ 

<T|| M -^ r || ioo(M n x[0iT]) <£ + c((r v ^) 2 + (/iv / ^) 2 ), Vx G M n , £ G [0,T]. 



The integral 



-(A-(t-rj)) 2 /(2?or 2 ) p -|x-fcm| 2 /(I>ft 2 +4i'A) 



^ m(x ' t) = j ~ (vh 2 +\uxyi 2 dX (5i 



cannot be taken analytically, but it allows obviously an approximate tensor product 
approximation. Making the substitution 

A * 



l + e-5 

we derive the integral over K 

OO 

f f e ~(rj-t/(l+^)) 2 /(V r 2 ) e -\x-hm\ 2 /(Vh 2 +4vt/(l+e-i)) 

K j>m {x,t) = -J — - - e _^ n/2 cosh 2^/ 2 ) 



with exponentially decaying integrand. Performing the last 2 substitutions in (3.5) 
we again transform the integrand to a doubly exponentially decaying function. 



16 



Approximations which converge with higher order to the solution of (5.1 ) can be 
obtained using quasi-interpolation of the right-hand side u by 



<A^t) = ^==^ E u(hm,rj) ms ( t -^-) nw(^) (5.9) 



Xi — hrrii " 



where t] 2 m are defined by (3.2). Since for all j < 2S — 1, ctj < 2M — 1 we have 

r(«+l)/2, j= ai = ...a n = 0, 



1=1 



TP 







otherwise . 



the following result can be derived in a standard way. 

Theorem 5.1. Given e > there exist T> > and T> > sncn £/iat for any 



r L fun 

oo 

estimate 



u G W^(M. n x M), with L = m&x(2M,2S), the quasi interpolant (5.9) satisfies the 



\ U (x,t) - ^ lT (x,t)| < Cl (m 2 ) M + c 2 (xv 2 ) 5 

t 2 ^ 1 (u /rVJal 2S - 1 
2_> ~\ \\ d x U \\L^(R^x[0,T}) + 
|a|=0 ' /3=0 



(5.10) 



'rJV n 



0\ 



t u \ |L°°(K n x[0,T]) 



where the constants c\ and c 2 do not depend on h, r, V, V and f . 
To obtain the cubature we use that 

t 



V2S 



T)2M 



X i 



s '-l) k {V Q r 2 ) k d 2k 



£ 

fc=0 
M 

£ 

fc=0 



k\ A k dt 2k 



-l) k (Vh 2 ) k d 2k ^ /(h2v) 



k\ A k dx 2 



2k 



Hence the heat potential of the quasi-interpolant (5.9) 



/ h)T (x,t) = / (P t _ A u fciT (.,A))(x)dA 



M 



7r(n+i)/2^p- l 
v u je 



u(hm, rj) ^fm ( x > *) 



k=0 



where we use the notation 

;-i)*(»„t 2 )* a 



fc=0 



s ■ ' ,Nfc (P 7- 2 )* <9 2fc n M ' "l \'>' / r> Z. - \ k ;)-/> 



i=l fc=0 



-l) fc (£>/i 2 ) fc ^ 



(5.11) 



k\ 4 fc 9a; 
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We have 

Q2k 



dx 



2k 



e -(* i -ftm i ) a /(©fc a -Hi/A) = ( l) fc /c!l fc _ ^ (-1/2)/^ - m*. 



(£>/i 2 + 4i/A)» 



£>/i 2 + 4z/A 



which by using (5.8) leads to the representation 



s / -i) k (v r 2 ) k d 2k 



k=0 



k\ A k dt 2k 



r (A-(^))V(^)-Q^ (A) ^_ mi)dA) 







admitting again a tensor product approximation. Here we denote by gu the function 

M 



g M {\x) = ^ 



2\fc 



(-1/2)/ a; 



k=0 



{Vh 2 + 4^A) fc+1 /2 fc V£>/i 2 + 4z/A 



-x 2 /(m 2 +4i/A) 



From Theorem 15. II it is easy to deduce that fh T approximates the solution / with 



the order 0{{\fD^T) 2S + (y/Vh) 2M ) plus the saturation error. 

Theorem 5.2. For any e > there exist T> > and T>q > Osuch that, for all 
u G W^(M. n x R), wrat/i L = max(2M, 25 1 ), i/je quasi-interp olant ( 5.11 ) approximates 
the solution of the Cauchy problem for the heat equation (5.1)-(5.2) with the error 
estimate 



|/(x,t) - / hlT (x,t)| < c 1>T (Vh 2 ) M + C 2 , T (V T 



.2\S 



2M-1 



(WZ>) 



l«l 



2S-1 



II«9"w||l°°(IR™x[0,T']) + ^ ^, ° ||df^|U°°(»"x[0,T]) 



ft! 



|a|=Q 13=0 

The constants c^t, i = 1,2, depend only on M and S. 
Proof. Since 



0\ 



1 



{A7Tu{t - X)) n / 2 

we obtain that 

|/(x,i) -f h)T {x,t)\ 



-|x-«|V(M*-a)) 



d£ = 1 



1 



dA 



-|x-€| a /(4i/(t-A)) 



(47n/) n / 2 y (t - a)™/ 2 

< T ||w - «fe, r ||i<x> (R n x[0! T]) , Vx G lR n , i G [0,T] 



|/(£,A)-A, T (£,A)|d£ 



From (5.10) the proof is complete. 



□ 
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5.1 Numerical example 

We have tested the approximation formula ( 5.5[ ) for solving the Cauchy problem 

ft~f xx = x 2 + t\ /(x,0) = 0, xeR, t>0, (5.12) 



t 3 

having the solution f(x,t)=t 2 -\ h tx 2 . 

3 



In the table [7] the difference fh jT (x,t) — f(x,t) for different values of h and r, 
T>q = T> = 2, at the time i = 0.01 and the point x = 0.01 is given. The numerical 
results in the table confirm that the error is 0(h 2 + r 2 ). 



r-^h- 1 


4 


8 


16 


32 


64 


128 


4 


1.25 10" 3 


7.81 10" 4 


6.64 10" 4 


6.35 10" 4 


6.27 10" 4 


6.25 10" 4 


8 


7.81 10~ 4 


3.12 10" 4 


1.95 10- 4 


1.66 10" 4 


1.58 10- 4 


1.56 10" 4 


16 


6.64 10" 4 


1.95 10- 4 


7.81 10- 5 


4.88 10" 5 


4.15 10" 5 


3.96 10" 5 


32 


6.34 10" 4 


1.66 HT 4 


4.88 10" 5 


1.95 10" 5 


1.22 10" 5 


1.03 10" 5 


64 


6.27 10" 4 


1.58 HT 4 


4.15 HT 5 


1.22 10" 5 


4.88 10" 6 


3.05 10- 6 


128 


6.25 10" 4 


1.56 10" 4 


3.96 10" 5 


1.03 10" 5 


3.05 10" 6 


1.22 10" 6 



Table 7: Error table for solving (5.12) with (5.5) 
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